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Abstract 

We examine the medium time quantum dynamics and population equilibration of two, three and 
four-well Bose-Hubbard models using stochastic integration in the truncated Wigner phase-space 
representation. We find that all three systems will enter at least a temporary state of equilibrium, 
with the details depending on both the classical initial conditions and the initial quantum statistics. 
We find that classical integrability is not necessarily a good guide as to whether equilibration will 
occur. We construct an effective single-particle reduced density matrix for each of the systems, 
using the expectation values of operator moments, and use this to calculate an effective entropy. 
Knowing the expected maximum values of this entropy for each system, we are able to quantify 
the different approaches to equilibrium. 

PACS numbers: 03.75.Lm, 03.75.Kk, 67.25.du, 03.75.Gg 



1 



I. INTRODUCTION 



The relaxation to equilibrium of closed quantum systems and their temporal evolution 



are important areas of study, as seen in, for example 



l|-|3|], with a beautiful experiment 



by Kinoshita et al. {4] having shown that relaxation to equilibrium does not happen in a 
trapped one dimensional Bose gas with contact interactions. This was not unexpected for 
a one dimensional untrapped Bose gas with point interactions, which is known to be an 
integrable system, but it had been thought that practical features such as the harmonic trap 
and imperfectly point-like interactions would compromise the integrability and the system 
would relax. On the other hand, there are closed quantum systems which are known to 
relax to a generalised equilibrium, without any interactions with a thermal cloud or other 
reservoir j2|, [sl, Isj. To the best of our knowledge, there is as yet no established consensus on 
the mechanism by which this happens. 

For the relaxation of closed quantum systems to a generalised equilibrium, one proposal 
is the eigenstate thermalisation hypothesis (ETH), in which every eigenstate of the Hamil- 
tonian implicitly contains a thermal state s], [t], with the initial coherence between the 
eigenstates being lost by dephasing during the dynamics. Srednicki, when introducing this 
hypothesis, claimed that a necessary condition was the validity of Berry's conjecture that 
the energy eigenfunctions behave like Gaussian random variables which is expected to 
hold for systems which exhibit classical chaos in at least a large majority of the classical 
phase space. 

In this work we study the relaxation or lack thereof in bosonic ultra-cold atomic systems 



held in two l9i , three 



10|, and four- well [Ul, ll2| tunnel-coupled potentials. For these sytems. 



the only constants of the motion are the total number of atoms and the total energy, so that in 
the classical sense we would expect only the twin well model to be integrable. Following the 
general wisdom, we would therefore not expect the two well system to relax to equilibrium, 
whereas we have previously shown that the four-well system, at least in one particular 
configuration, will demonstrate this feature over medium times 12| . 

We use stochastic integration in the truncated Wigner representation [isl, [l^ to perform 
phase-space analyses of the quantum dynamics and calculate the expectation values of the 
well populations and other operator products which allow us to construct effective reduced 
single-particle density matrices. From these matrices we are able to calculate a pseudoen- 
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12 1 . For the two-site model 



tropy which we previously found useful for the four well system 
we are also able to calculate exact quantum results using a matrix equation for the number 
state coefficients, finding excellent agreement with the stochastic results over most of the 
time domain. We also calculate a type of Lyupanov exponent classically, using coupled 
Gross-Pitaevskii equations, and investigate the extent to which this allows us to predict 
which systems will relax. We find that all the systems we investigate will reach a gener- 
alised, but not necessarily static, equilibrium, but not for arbitrary initial quantum states. 
We find that the differences can be predicted qualitatively via the off-diagonal elements of 
the reduced density matrices at the beginning of the time evolution. 



II. PHYSICAL MODEL, HAMILTONIAN AND EQUATIONS OF MOTION 

We consider three different physical arrangements, firstly with two wells of equal depth, 
secondly with three equal wells at the vertices of an equilateral triangle, and thirdly with 
four equal wells in a square arrangement. We will follow a procedure which is effectively 
equivalent to changing the Hamiltonians at t = 0, in that all our initial conditions will have 
at least one well unoccupied. We outline here the the standard procedure for developing an 



effective Hamiltonian for two wells {q], where we consider separate wells with an independent 
condensate in each of the at the beginning of our investigations. The procedures for three 



or four wells are a simple extension of this [12|. The Hamiltonian for a condensate in an 
external trapping potential, Vext{f^, may be written as 



n= dr 



2m 



(1) 



where ip is the field operator for the condensate, and the non-linear interaction parameter 
is Uq = 2Tcah^ /m, where a is the s-wave scattering length describing two-body collisions 
within the condensate, and m is the atomic mass. In the case where the external potential 
provides a two well confinement for the condensate, we may simplify the above Hamiltonian 
by making use of the two-mode approximation. At zero temperature all atoms in the system 
are condensed and if the ground state energies of the condensate in either well are sufficiently 
separated from the energies of the condensate in all other excited single particle states, 
transitions to or from the modes of interest and these higher lying states can be neglected. 
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We may then expand the field operator as 



(2) 



where the dj are bosonic annihilation operators in each of the wells, and the 0j are the 
ground state spatial wave functions of the condensate in each of the wells. 
Using this in Eq. ([1]), we find an effective Hamiltonian 

^ ^ -j- ^ ^ -j- ^ 

Heff = EihiCii + E2(l2<i2 

+hx ^d|d|didi + a\d\a2a2^ 

—hJ (^a[a2 + alcii^ , (3) 

where we have neglected the spatial overlap of the different well densities. The single well 
bound state energies, Ei, are 

Ei = J df<p*{r) (^V^ + VU^^ Ur)- (4) 
J, the tunnel coupling, is 

J=^ j dr<Plif^ {^V' + K.*(r)) 02(r1, (5) 
and effective non-linear interaction term is 

X = df\Ur)t. (6) 

We set the single well bound state energies equal because we will consider only symmetric 
potentials where we can set Ei = E2 = 0. 

This process simplifies the analyses by approximating the atoms in each separate well as 
being in a single mode, meaning that our equations are much easier to solve than would 
otherwise be the case. Generalising to the three and four site models, we will use {i = 
1,2,3) as bosonic annihilation operators in each of the two or three wells, and Sj and bi, 
with i = 1,2, in each side of the four-well system. We parametrise time by setting J = 1, 
so that dimensionless time as displayed in the results is labelled as Jt. A schematic of the 
three well configuration is shown in Fig. [H with the four well system being treated as shown 
in Fig. H 
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FIG. 1: Schematic of our three- mode Bose-Hubbard system. The hi are the bosonic annihilation 
operators for each mode, while J represents the coupling rate between the modes. In this article, 
we always set J = 1, which sets the units of time. 



We may now also write the effective Hamiltonians for the three and four site systems, 
generalising from that for twin wells. For the three-well system we find the effective Hamil- 
tonian, 

3 
i=l 

—hJ {^\a2 + agO-i + ajds + a\ai + dgO-a + 030-2 j . (7) 

Finally, for the the four-well system we have, 
2 

1=1 

—hJ {a\a2 + 02^1 + h\h2 + h^bi + a\hi + h\ai + a\h2 + ^2^2^ • (8) 

In order to calculate the quantum d yna mics of these three systems, we will make use 
of the truncated Wigner representation jisl. \\\. which gives results indistinguishable from 
those of the full quantum matrix equations for the number state coefficients p^] up until at 
least the relaxation time in the case of two and three wells, and is much easier to use for four- 
wells, where the full matrix equations become very difficult to use. Following the standard 

n 

methods pj( , we find generalised Fokker-Planck equations for the Wigner pseudoprobability 
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FIG. 2: Schematic of our four-mode Bose-Hubbard system. The en and hi are the bosonic annihi- 
lation operators for each mode, while J represents the coupling rate between the modes. 

functions. As these contain third order derivatives, they cannot be mapped onto stochastic 
differential equations. Although methods have been developed to allow a mapping onto 
stochastic difference equations the numerical integration of these equations is even 



more unstable than that of the positive-P representation equations [Ij, Il8|], so we will 
not follow this route here. Instead, we truncate the third order terms in the generalised 
Fokker-Planck equations and make a mapping to coupled equations for the Wigner variables 
corresponding to the operators found in the Hamiltonians. We note here that classical 
averages of the Wigner variables correspond to symmetrically ordered operator expectation 
values, so that the necessary reordering must be undertaken before we arrive at solutions 
for physical quantities, for which normal ordering is more appropriate. 

Having done this, we find the sets of coupled equations given immediately below. For the 
two-well system, the truncated Wigner equations of motion are, 

dai 



dt 
da2 



—2ix\oii\'^ ai + iJa2 



(9) 
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for the three-well system they are 

dai 

~dt 

da2 

~dt 

da^ 

~~dF 

and for the four-wells they are 

dai 

~dt 

da2 

H 

dl3i _ 

~dt 

dA ^ 
dt 



-2zx|ai|^ai + iJ (02 -l- as) 
-2ix\a2\^a2 + iJ (ai + as) 
-2zx|as|^a3 + iJ (ai + a2) 



-2ix\ai\^ai + iJ {a2 + A) 
-2ix|a2pa2 + iJ (ai + 132) 
-21x^131 + 13 {132 + 
-2ix\f32?f32 + iJ{l3i + a2). 



(10) 



(11) 



Although the three sets of equations above might look classical, the Wigner variables 
themselves are drawn from appropriate distributions for the desired initial quantum states, 
with the stochasticity coming from the initial conditions for the chosen quantum states. The 
truncated Wigner equations are solved numerically by taking averages over a 
of stochastic trajectories, with initial conditions sampled probabilistically 
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arge number 



19|. 



III. CLASSICAL INTEGRABILITY AND STABILITY 

Before we investigate the quantum dynamics further we will examine the classical systems, 
which are described by coupled Gross-Pitaevskii equations which have the same form as 
above, but have deterministic initial conditions. Berry's conjecture [8] states that a quantum 
system can thermalise if the related classical system is unstable or chaotic in at least a 
large majority of the classical phase space. To examine the applicability of this conjecture, 
we numerically calculate approximate Lyupanov exponents for the three and four mode 
systems [201], which we define for the dynamics of one of the coupled wells as 

_ j.^ l M5N,{r)) 



where 



SN^{T) = \M'\T)~Nf{r)l J = 1,2,3,4, (13) 
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where A^j^^ is an initial condition slightly perturbed from Nj^\ Note that we have used 
numbers here rather than complex amplitudes, as it is the numbers which are directly 
observable. In practice, we obviously cannot integrate the equations for infinite time, so 
we integrate the coupled GPE type equations over a reasonably long time and look at 
the development of SNj{t) and hence Lj{t). We also make the caveat here that because 
the system is both bounded and periodic, these effective Lyupanov exponents are not as 
definitive as in a non-periodic system. 

In these systems as we define them, the only constants of the motion are the total energy 
and the total number of atoms. The obvious classical degrees of freedom are the number 
of atoms in each well, which, given the constraint on total number, would seem to be one 
less than the number of wells. This approach suggests that the three-well system has the 
same number of classical degrees of freedom as it does constants of the motion, and therefore 
should be classically integrable. The four-well system, as we have previously shown, is not 
stable, at least when two different tunneling rates are present 12|. 

In Fig. I ^a) we show the exponents for the three-well system, beginning in a situation 
far from the expected equilibrium, which would have one third of the atoms in each well. 
The classical solutions are again regular and periodic, with full oscillations between and 
300 atoms in well three, with the other two oscillating between 150 and 0. The perturbed 
classical solutions, while still regular and periodic, never see more than 149 or less than 1 
atom in wells one and two, while well three oscillates between 2 and 302. Fig. [ ^b)| shows 
that the four-well system also behaves in a similar manner. 



IV. RELAXATION 



In a state of relaxed equilibrium, we expect that the entropy of the system will be max- 
imised. In a quantum system, the entropy is normally defined using the density matrix for 
the system. While this works very well for small systems, such as those that are of interest 
in discrete variable quantum information applications, it becomes difficult to calculate the 
full density matrix for a many-body interacting quantum system, which will have a much 
larger Hilbert space. As in our previous work [12j, we can define something which behaves 
as a reduced single-particle density matrix, which then allows us to calculate an effective 
entropy. Importantly for possible experimental investigations, all the quantities needed are 
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FIG. 3: (colour online) The time development of the Lyupanov exponents for the three- and four- 



well systems. For (a) the reference initial condition was with all 300 atoms equally distributed 
between wells one and two, with none in well three, and the perturbed initial condition had 149 
in wells one and two, with 2 atoms in well three, x^tot = 0.5 J, with J = 1. The solid line is 
Li, the dash-dotted is L2 (dynamics indistinguishable from Li) and the dotted line is L3. For |(b)| 
the reference initial condition was with 200 atoms in each of wells ai and 621 with the other two 
unoccupied. The perturbed initial condition had 199 in wells oi and b2, with 1 atom in each of 
the other two. x^tot = 0.5 J, with J = 1. The four exponents are indistinguishable. All quantities 
plotted in this and subsequent plots are dimensionless. 

rn 

in principle measurable using techniques developed by Ferris et al. |21j]. For the four- well 
system we define, 

^ {a[ai) {a[a2) {a[hi) {a[h2) \ 

1 (0.20-1) (02^2) 

Pa = — 

{Ntot) {h[ai) {})%) {blh) {b\b2) 

\{bla^) {bla2) {blh) (b\b2) J 

with the matrices for two and three wells being obvious truncations of the above. It is then 

an easy matter to calculate a single-particle pseudoentropy from this matrix. 



(14) 



e,(t) = -Tr[p,(t) lnp,(t)]. 



(15) 



which will have a maximum value of In j when the atoms are equally distributed throughout 
the j wells, which is statistically the most probable situation. 
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FIG. 4: (colour online) The mean-field population dynamics and pseudoentropy for the two-well 
system, with A^i(O) = 200, A''2(0) = 0, x^tot = 0.5 J, and J = 1 , calculated using the truncated 
Wigner representation. The solid (blue) lines are the average over 1.57 x 10^ trajectories for an 
initial Fock state, while the dash-dotted lines (red) are averaged over 2.33 x 10^ trajectories for an 
initial coherent state, (a) gives the expectation values of the population in well one and [(b)] gives 



the calculated pseudentropies. 

In Fig. Hlwe show the results for the two- well system, beginning with all 200 atoms initially 
in one of the wells. From Fig. [ ^a)| we see that the atoms equalise between the two wells, 
with the initial quantum states having almost no discernible effect on the process. This is 
also clear when we look at Fig. I ^b) where the pseudoentropy rises to a value very close to 



In 2 ^ 0.6931, which is the maximum single-particle value expected of this system in thermal 
equiUbrium. We therefore see that, once quantum effects are included, this system appears 
to relax to something close to its equilibrium state at zero temperature without any contact 
with a reservoir, despite the fact that it is classically integrable. We have also integrated the 
exact quantum equations for the number state coefficients using a matrix method for initial 



Fock states 
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22| and see a partial revival of the oscillations centred around Jt = 420, as 



shown in Fig. [^a)[ We used this method out to Jt = 900, finding one more very partial 



revival at around Jt = 800. In Fig. [ ^b) we show the P{Ni) at Jt = 900 for well one of 
this system. The true equilibrium ground state for repulsive interatomic interactions would 
be a binomial distribution between the two wells, with the distribution in each well closely 
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FIG. 5: (colour online) The first mean-field partial revival and the probablities for each number 
in well one at Jt = 900, calculated using the exact method for the same parameters and initial 
conditions as in Fig. HI but with initial Fock states. 

approximating a Gaussian centred on Nj = 100. We see that this is not the case here, with 
the highest probably actually being for all the atoms in well one. This shows that, despite 
the pseudoentropy being close to its maximum and the mean populations being equal at 
this time, the system is not in a true equilibrium state. We will return to this issue below. 

When we look at the triple well system, we find differences which depend entirely on the 
quantum statistics, in a similar manner to previous investigations of quantum superchem- 



istry 
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24( 1 ■ As shown in Fig. |6l when we start with a total of 300 atoms, half in one well 
and half in another, with the third well unoccupied, the subsequent behaviours for initial 
coherent and Fock states are qualitatively different. The initial Fock states equilibriate to 
an equal number in each well and the pseudoentropy climbs to very close to In 3 ~ 1.0986, 
but initial coherent states evolve to a different final distribution. We find approximately 87 
atoms in each the wells which were originally occupied, while the initially unoccupied well 
holds a final number of 126. The final value of the pseudoentropy is well below the possible 
maximum, as the state of the three wells at the end of our integration time is by no means 
the statistically most probable distribution. A clue to the explanation of this behaviour can 
be found in the reduced single-particle density matrices at the final time of Jt = 200. For 
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200 
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FIG. 6: (colour online) The mean-field population dynamics of well 1 and the pseudoentropy for 
the same parameters and initial conditions as the reference state in Fig. [ ^a)| in the three-well 
model. The solid (blue) lines are the average over 1.04 x 10^ trajectories for initial coherent states, 
while the dash-dotted lines (red) are averaged over 8 x 10^ trajectories for initial Fock states. The 
total atom number was 300, with half of these initially in well 1 and half in well 2. 



gives the 



expectation values of the population in well one and |(b)| gives the calculated pseudoentropies. 



the initial Fock states we find 
/ 



P3 



v 



0.3315 0.0031 - 0.0017i 

0.0031 + 0.0017Z 0.3297 
-0.0023 - O.OOOli -0.0019 + 0.0019i 



-0.0023 + 0.0001i 
-0.0019 - 0.0019i 
0.3388 



(16) 



while for initial coherent states the final result is 
/ 



P3 



0.2895 0.2880 + O.OOOOi 0.1039 + 0.0004i 
0.2880 - O.OOOOi 0.2897 0.1039 + 0.0003i 
0.1039 - 0.0004i 0.1039 - 0.0003i 0.4209 



(17) 



from which we see that the off-diagonal elements are much larger for the initial coherent 
states. As these elements represent coherences between the modes, they contain information 
and therefore, the larger they are, the further the system is from the maximum entropy 
equiUbrium state. 

For the four-well system, as shown in Fig. [TJ we see that the maximum pseudoentropy 
(In 4 ^ 1.3863) is closely approached for both initial Fock and coherent states, but also that 
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150 



FIG. 7: (colour online) The population dynamics of well aiof the four-well system and the pseu- 
doentropy for the same parameters and initial conditions as the reference state of Fig. [ ^b)[ The 
solid (blue) lines are the average over 1.36 x 10^ trajectories for initial Fock states, while the dash- 
dotted lines (red) are averaged over 4.05 x 10^ trajectories for initial coherent states. The total 



atom number was 400, with half of these initially in well ai and half in well 62- (a) gives the 



expectation values of the population in well ai and |(b)| gives the calculated pseudoentropies. 

the population dynamics are quite different. The small oscillations about the equilibrium 
value in well oi die down much slower for initial coherent states than for initial Fock states. 
The off-diagonal elements of the matrix at the end of the integration time were also larger 
for initial coherent states, although the number distribution is obviously tending towards a 
quarter of the total in each well. 

We note here that for initial Fock states in all these systems, the off-diagonal elements of 
the reduced density matrix are all zero at t = 0, while this is not the case for initial coherent 
states. This can be easy seen by considering some of the typical initial matrix elements for 
each system. For the two well system with an initial coherent state in one side and nothing 
in the other, we find 

(«!, 0|a|a2|«i, 0) = {ai,0\d\ai\ai,0) = 0, (18) 

while for a Fock state in one side and nothing in the other, we have 

{Ni,0\ald2\Ni,0) = (Ari,0|a^ai|Ari,0) = 0, (19) 

so that all the off-diagonal elements are originally zero. For the three-well system with our 
coherent state initial conditions, we find two of the off-diagonal elements are non-zero, 

(cti, a2, 0|ajd2|ai, 02, 0) = ala2 = {{cti, 02, 0|a{a2|Q;i, a2, 0))*, (20) 
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with expectation values of any operator products containing or al being zero. For the 
initial Fock states in wells one and two, we find that all the off-diagonal elements are again 
initially zero. With the four-well system and the two Fock states, all off-diagonal elements 
are initially zero, whereas for the two coherent state occupations, we find the non-zero 
elements, 

{ai,0,0,/32\a[b2\ai,0,0,/32) =al/32 = ((ai, 0, 0, /32|6^ai|ai, 0, 0, ^Ss))*- (21) 

We therefore see that, where the initial reduced density matrices are the same for each choice 
of quantum states, as in the two-well example, we see almost indistinguishable evolution of 
both the populations and the pseudoentropy. Where they exhibit the most difference, as 
in the three-well example, the subsequent dynamics are very different. The four-well case, 
where the initial matrices are less different, lies between these two extremes. 

V. CONCLUSIONS AND DISCUSSION 

Our stochastic phase-space investigations of the quantum dynamics of these three different 
Bose-Hubbard models shows that all three may equilibrate for given initial conditions. We 
have also shown that an effective reduced single-particle density matrix may be calculated 
using phase-space methods. While this matrix does not have all the information contained in 
the full quantum density matrices, it does allow for the calculation of a pseudoentropy and an 
investigation of the equilibriation process. Perhaps more importantly, all the measurements 
necessary to construct this density matrix are possible in principle. We have shown that 
classical integrability is not necessarily a good guide as to whether a given system will relax or 
not, as seen by the two- well system, which does relax to a close to equilibrium state, with only 
minor revivals over the time we investigated. We have also shown that the presence of initial 
coherences between the different modes can affect the equilibriation process, as shown by 
the presence of the off-diagonal elements of the density matrix. We have also demonstrated 
that the stochastic phase-space methods developed for quantum optics can be useful for the 
theoretical study of statistical mechanical processes in interacting atomic systems where the 
Hilbert space becomes so large as to make density matrix methods extremely difficult to 
use. 

What we cannot show, due to the limitations of our numerical methods, is whether the 
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oscillations in these systems ever revive fully in finite time, after their initial collapse. We 
suffer from two restrictions here. The first is that our exact method, although applicable 
to the two and three well models, can still not be used for arbitrary times. The second 
is that the truncated Wigner method gives us accurate information on the collapse of the 
oscillations, but not on the revivals. However, taken together, these two methods allow 
us to investigate the medium time dynamics and also allow us to determine the effects of 
different initial quantum states on the mean-field dynamics. In none of the systems we have 
examined here do the solutions of the mean-field classical equations even approximate the 
true dynamics of the mean-fields. In that sense these are truly quantum systems. On a final 
note, it is interesting to consider the systems we have examined as having non-Markovian 
reservoirs connected to them at t = 0, when the unoccupied wells become accessible to the 
atoms. Whether the relaxation we see can be described accurately as a smaller subsytem 
coming into some type of equilibrium with this bath will be investigated in later work. 
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